GLMM Part2

Preparations

Load the necessary libraries

library(car)       #for regression diagnostics
library(broom)     #for tidy output
library(ggfortify) #for model diagnostics
library(sjPlot)    #for outputs
library(knitr)     #for kable
library(effects)   #for partial effects plots
library(ggeffects) #for partial effects plots
library(emmeans)   #for estimating marginal means
library(MASS)      #for glm.nb
library(MuMIn)     #for AICc
library(tidyverse) #for data wrangling
library(nlme)
library(lme4)
library(glmmTMB)
library(broom.mixed)
library(glmmTMB)   #for glmmTMB
library(DHARMa)   #for residuals and diagnostics
library(performance) #for diagnostic plots
library(see)         #for diagnostic plots

Scenario

To investigate differential metabolic plasticity in barramundi (Lates calcarifer), Norin, Malte, and Clark (2015) exposed juvenile barramundi to various environmental changes (increased temperature, decreased salinity and increased hypoxia) as well as control conditions. Metabolic plasticity was calculated as the percentage difference in standard metabolic rate between the various treatment conditions and the standard metabolic rate under control conditions. They were interested in whether there was a relationship between metabolic plasticity and typical (control) metabolism and how the different treatment conditions impact on this relationship.

A total of 60 barramundi juveniles were subject to each of the three conditions (high temperature, low salinity and hypoxia) in addition to control conditions. Fish mass was also recorded as a covariate as this is known to influence metabolic parameters.

Barramundi

Sampling design

Format of norin.csv data files

FISHID	MASS	TRIAL	SMR_contr	CHANGE
1	35.69	LowSalinity	5.85	-31.92
2	33.84	LowSalinity	6.53	2.52
3	37.78	LowSalinity	5.66	-6.28
..	..	..	..	..
1	36.80	HighTemperature	5.85	18.32
2	34.98	HighTemperature	6.53	19.06
3	38.38	HighTemperature	5.66	19.03
..	..	..	..	..
1	45.06	Hypoxia	5.85	-18.61
2	43.51	Hypoxia	6.53	-5.37
3	45.11	Hypoxia	5.66	-13.95

FISHID	Categorical listing of the individual fish that are repeatedly sampled
MASS	Mass (g) of barramundi. Covariate in analysis
TRIAL	Categorical listing of the trial (LowSalinity: 10ppt salinity; HighTemperature: 35 degrees; Hypoxia: 45% air-sat. oxygen.
SMR_contr	Standard metabolic rate (mg/h/39.4 g of fish) under control trial conditions (35 ppt salinity, 29 degrees, normoxia)
CHANGE	Percentage difference in Standard metabolic rate (mg/h/39.4 g of fish) between Trial conditions and control adjusted for 'regression to the mean'.

Read in the data

norin = read_csv('../public/data/norin.csv', trim_ws=TRUE)
glimpse(norin)

## Rows: 180
## Columns: 5
## $ FISHID    <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
## $ MASS      <dbl> 35.69, 33.84, 37.78, 26.58, 37.62, 37.68, 30.62, 50.37, 24.…
## $ TRIAL     <chr> "LowSalinity", "LowSalinity", "LowSalinity", "LowSalinity",…
## $ SMR_contr <dbl> 5.847466, 6.530707, 5.659556, 6.278200, 4.407336, 4.818589,…
## $ CHANGE    <dbl> -31.919389, 2.520929, -6.284968, -4.346675, -3.071329, -15.…

Since the response variable is a difference between the SMR under control conditions compared to the specific trial conditions (either Low salinity, High temporature or Hypoxia), the only possible distribution we can use is a Gaussian. Well to clarify this, we could use a t-distribution, but this is not really available for frequentist mixed effects routines in R.

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of temperature and (centered) mean fish size on SDA peak. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual fish.

Lets begin by declaring the categorical predictors and random effects as factors.

norin = norin %>% mutate(FISHID=factor(FISHID),
                         TRIAL=factor(TRIAL))

Each of the fish were exposed to each of the three trials (treatments). So lets explore the distributions of responses within each of these trials.

ggplot(norin, aes(y=CHANGE, x=TRIAL)) + geom_boxplot()

Conclusions:

• these seem reasonable enough

The researchers considered that physiological plasticity might be effected by the fishes basal metabolic rate. For example, a fish with relatively high metabolism might have less scope to increase this metabolism than a fish with a lower metabolism. Therefore, the SMR under control conditions (just prior to the specific trial) could be added as a continuous covariate.

Doing so would introduce two additional assumptions:

1. linearity: that the relationship between the response and control SMR is linear (if we intend to model it as a linear trend)
2. homogeneity of slopes: that the rate of change in response associated with a one unit change in control SMR is similar for each trial (unless we include an interaction)

ggplot(norin, aes(y=CHANGE, x=SMR_contr, shape=TRIAL, color=TRIAL)) +
    geom_smooth(method='lm') + geom_point()

ggplot(norin, aes(y=CHANGE, x=SMR_contr, shape=TRIAL, color=TRIAL)) +
  geom_smooth() + geom_point()

It might also be worth exploring how consistent the trial effect is across the different fish. This can give us an idea of whether the addition of a random intercept/slope model would be useful.

ggplot(norin, aes(y=CHANGE, x=as.numeric(FISHID), color=TRIAL)) +
    geom_point() + geom_line()

Conclusions:

• generally, the difference between control and trial SMR is greatest for the High temperature trial.
• there is less consistency about the relative effects of Hypoxia and Low Salinity.
• hence a random intercept/slope model might be useful.

Finally, as an acknowledgement that the metabolic response might be influenced

by the mass of the fish, the researchers contemplated including fish Mass as a continuous covariate. There are numerous alternative ways to incorporate covariates that are known to impact on a response.

• add them as an offset. This will include the covariate in the model, yet will not estimate a partial slope for this predictor. Rather the slope is assumed to be 1. Since the parameter does not need to be estimated, adding a predictor as an offset does not consume any degrees of freedom. However, it does assume that there is a slope of one. This approach can be useful as an alternative to using the response divided by the covariate as the response. That is, it is equivalent to using a response that has been standardized by the covariate, yet still allowing the modelling distribution that would be appropriate for the original response.
• add them as a regular covariate. This will include the covariate in the model and esimtate a partial slope parameter just like any other term. Although it will incur a degrees of freedom penalty, it does not assume a slope of 1 and could be modelled as non-linear if appropriate.

Lets explore the relationship between the response and Mass separately for each Trial.

ggplot(norin, aes(y=CHANGE, x=MASS, color=TRIAL)) +
  geom_point() +
  geom_smooth(method='lm')

Conclusions:

• there does not seem to be much of a relationship between the response and Mass. Nevertheless, in the presence of control SMN, it might still be a useful covariate at explaining some of the unexplained variability (and thus increasing the power of the model).
• there is no evidence that the response and Mass are related in a non-linear manner.

Fit the model

For each of the following modelling alternatives, we will:

• establish whether MASS is a useful covariate (based on AICc). This step involves evaluating alternative fixed effects structures. When doing so, we must use Maximum Likelihood (ML) rather than Residual Maximum Likelihood (REML).
• establish whether a random intercept/slope model is useful (based on AICc). This step involves evaluating alternative random effects structures. When doing so, we should use Restricted Maximum Likelihood (REML)

lme (nlme)

fixed structures

##Compare models that estimate partial slope for MASS vs an offset for MASS
##must use ML to compare models that vary in fixed effects
norin.lme1 <- lme(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + offset(MASS), random=~1|FISHID, data=norin, method='ML')
## Unfortunately,  update() does not remove offset().
## We will just have to write the other model out in full as well.
norin.lme2 <- lme(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + MASS,  random=~1|FISHID, data=norin, method='ML')
## Now without MASS altogether
norin.lme3 <- lme(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE),  random=~1|FISHID, data=norin, method='ML')

## Compare these models via AICc
AICc(norin.lme1,norin.lme2, norin.lme3)

## Alternatively,  we can use sequential Likelihood Ratio Tests (LRT)
anova(norin.lme1, norin.lme2, norin.lme3)

Conclusions:

• on the basis of AICc, the model without MASS is considered the most parsimonious (lowest AICc)
• we will proceed with the fixed structure of model 3

random structures

norin.lme3a = update(norin.lme3, method='REML')
norin.lme3b = update(norin.lme3a, random=~TRIAL|FISHID)
AICc(norin.lme3a,  norin.lme3b)

anova(norin.lme3a, norin.lme3b)

Conclusions:

• the AICc of the random intercept/slope model is lowest
• therefore we will proceed with a random intercept/slope model

lmer (lme4)

fixed structure

##Compare models that estimate partial slope for MASS vs an offset for MASS
##must use ML to compare models that vary in fixed effects
norin.lmer1 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + offset(MASS) + (1|FISHID), data=norin, REML=FALSE)
## Unfortunately,  update() does not remove offset().
## We will just have to write the other model out in full as well.
norin.lmer2 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + MASS + (1|FISHID), data=norin, REML=FALSE)
## Now without MASS altogether
norin.lmer3 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE)+ (1|FISHID), data=norin, REML=FALSE)

## Compare these models via AICc
AICc(norin.lmer1,norin.lmer2, norin.lmer3)

## Alternatively,  we can use sequential Likelihood Ratio Tests (LRT)
anova(norin.lmer1, norin.lmer2, norin.lmer3)

Conclusions:

• on the basis of AICc, the model without MASS is considered the most parsimonious (lowest AICc)
• we will proceed with the fixed structure of model 3

random structure

norin.lmer3a = update(norin.lmer3, REML=TRUE)
## unfortunately,  the following random intercept/slope model cannot be fit
##norin.lmer3b = update(norin.lmer3a, ~ . -(1|FISHID) +(TRIAL|FISHID))
##AICc(norin.lmer3a,  norin.lmer3b)
##anova(norin.lmer3a, norin.lmer3b)

Conclusions:

• unfortunately, we are unable to fit a model with random intercept/slope
• therefore we will proceed with a random intercept model

glmmTMB

fixed structure

##Compare models that estimate partial slope for MASS vs an offset for MASS
norin.glmmTMB1 = glmmTMB(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + offset(MASS) + (1|FISHID), data=norin, REML=FALSE)
## Unfortunately,  update() does not remove offset().
## We will just have to write the other model out in full as well.
norin.glmmTMB2 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + MASS + (1|FISHID), data=norin, REML=FALSE)
## Now without MASS altogether
norin.glmmTMB3 <- lmer(CHANGE ~ TRIAL*scale(SMR_contr, scale=FALSE) + (1|FISHID), data=norin, REML=FALSE)

## Compare these models via AICc
AICc(norin.glmmTMB1,norin.glmmTMB2, norin.glmmTMB3)

Conclusions:

• on the basis of AICc, the model without MASS is considered the most parsimonious (lowest AICc)
• we will proceed with the fixed structure of model 3

random structure

norin.glmmTMB3a <- update(norin.glmmTMB3, REML=TRUE)
## unfortunately,  the following random intercept/slope model cannot be fit
## norin.glmmTMB3b <- update(norin.glmmTMB3a, ~ . - (1|FISHID)+ (TRIAL|FISHID))
##AICc(norin.glmmTMB1a,  norin.glmmTMB1b)

Conclusions:

• unfortunately, we are unable to fit a model with random intercept/slope
• therefore we will proceed with a random intercept model

Model validation

lme (nlme)

Conclusions:

• whilst the residual plot is ok, the Q-Q plot deviates substantially from a straight line. This suggests that there is more mass in the tails of the distribution than a Gaussian would expect. Although the model is reasonably robust to this, it would be nice to have the option to model these data against a flatter distribution such as a Student t. Unfortunately, this is not available.

plot_model

plot_grid(plot_model(norin.lme3b, type='diag')[-2])

Conclusions:

• whilst the residual plot is ok, the Q-Q plot deviates substantially from a straight line. This suggests that there is more mass in the tails of the distribution than a Gaussian would expect. Although the model is reasonably robust to this, it would be nice to have the option to model these data against a flatter distribution such as a Student t. Unfortunately, this is not available.

Performace model checking

performance::check_model(norin.lme3b)

## Could not compute standard errors from random effects for diagnostic plot.
## Homogeneity of variance could not be computed. Cannot extract residual variance from objects of class 'lme'.

DHARMa residuals

#norin.resid = simulateResiduals(norin.lme3a,  plot=TRUE)

lmer (lme4)

plot_model

plot_grid(plot_model(norin.lmer3a, type='diag')[-2])

Conclusions:

• whilst the residual plot is ok, the Q-Q plot deviates substantially from a straight line. This suggests that there is more mass in the tails of the distribution than a Gaussian would expect. Although the model is reasonably robust to this, it would be nice to have the option to model these data against a flatter distribution such as a Student t. Unfortunately, this is not available.

Performace model checking

performance::check_model(norin.lmer3a)

Conclusions:

• whilst the residual plot is ok, the Q-Q plot deviates substantially from a straight line. This suggests that there is more mass in the tails of the distribution than a Gaussian would expect. Although the model is reasonably robust to this, it would be nice to have the option to model these data against a flatter distribution such as a Student t. Unfortunately, this is not available.

DHARMa residuals

norin.resid = simulateResiduals(norin.lmer3a,  plot=TRUE)

Conclusions:

• in contrast to the regular Q-Q normal plots, the Q-Q normal plots generated from the simulated residuals look fine.
• the residual plot from the simulated residuals look fine
• the diagnostic tests reveal no issues.

glmmTBM (glmmTMB)

plot_model

plot_grid(plot_model(norin.glmmTMB3a, type='diag')[-2])

Conclusions:

• whilst the residual plot is ok, the Q-Q plot deviates substantially from a straight line. This suggests that there is more mass in the tails of the distribution than a Gaussian would expect. Although the model is reasonably robust to this, it would be nice to have the option to model these data against a flatter distribution such as a Student t. Unfortunately, this is not available.

Performace model checking

performance::check_model(norin.glmmTMB3a)

Conclusions:

• whilst the residual plot is ok, the Q-Q plot deviates substantially from a straight line. This suggests that there is more mass in the tails of the distribution than a Gaussian would expect. Although the model is reasonably robust to this, it would be nice to have the option to model these data against a flatter distribution such as a Student t. Unfortunately, this is not available.

DHARMa residuals

norin.resid = simulateResiduals(norin.glmmTMB3a,  plot=TRUE)

Conclusions:

• in contrast to the regular Q-Q normal plots, the Q-Q normal plots generated from the simulated residuals look fine.
• the residual plot from the simulated residuals look fine
• the diagnostic tests reveal no issues.

Partial plots model

lme (nlme)

plot_model

plot_model(norin.lme3b,  type='eff',  terms=c('SMR_contr', 'TRIAL'),  show.data=TRUE)

We can also explore the fixed effects (as a caterpillar plot).

plot_model(norin.lme3a,  type='est')

And similarly, the individual random effects (intercepts of each fish).

#plot_model(norin.lme3a,  type='re')

allEffects

plot(allEffects(norin.lme3b), multiline=TRUE,  ci.style='bands')

ggpredict

ggpredict(norin.lme3b,  c('SMR_contr', 'TRIAL')) %>% plot(add.data=TRUE)

ggemmeans

#ggemmeans(norin.lme3a,  ~SMR_contr*TRIAL) %>% plot(add.data=TRUE)

lmer (lmer)

plot_model

plot_model(norin.lmer3a,  type='eff',  show.data=TRUE) %>% plot_grid

plot_model(norin.lmer3a,  type='eff',  terms=c('SMR_contr', 'TRIAL'),  show.data=TRUE)

We can also explore the fixed effects (as a caterpillar plot).

plot_model(norin.lmer3a,  type='est')

And similarly, the individual random effects (intercepts of each fish).

plot_model(norin.lmer3a,  type='re')

allEffects

plot(allEffects(norin.lmer3a), multiline=TRUE,  ci.style='bands')

ggpredict

ggpredict(norin.lmer3a,  c('SMR_contr', 'TRIAL')) %>% plot(add.data=TRUE)

ggemmeans

ggemmeans(norin.lmer3a,  ~SMR_contr*TRIAL) %>% plot(add.data=TRUE)

glmmTMB (glmmTMB)

plot_model

plot_model(norin.glmmTMB3a,  type='eff',  show.data=TRUE) %>% plot_grid

plot_model(norin.glmmTMB3a,  type='eff',  terms=c('SMR_contr', 'TRIAL'),  show.data=TRUE)

We can also explore the fixed effects (as a caterpillar plot).

plot_model(norin.glmmTMB3a,  type='est')

And similarly, the individual random effects (intercepts of each fish).

plot_model(norin.glmmTMB3a,  type='re')

allEffects

plot(allEffects(norin.glmmTMB3a), multiline=TRUE,  ci.style='bands')

ggpredict

ggpredict(norin.glmmTMB3a,  c('SMR_contr', 'TRIAL')) %>% plot(add.data=TRUE)

ggemmeans

ggemmeans(norin.glmmTMB3a,  ~SMR_contr*TRIAL) %>% plot(add.data=TRUE)

Model investigation / hypothesis testing

lme (nlme)

summary

summary(norin.lme3b)

## Linear mixed-effects model fit by REML
##  Data: norin 
##        AIC      BIC    logLik
##   1593.961 1635.028 -783.9803
## 
## Random effects:
##  Formula: ~TRIAL | FISHID
##  Structure: General positive-definite, Log-Cholesky parametrization
##                  StdDev    Corr         
## (Intercept)      16.323848 (Intr) TRIALH
## TRIALHypoxia     14.926864 -0.601       
## TRIALLowSalinity 30.304028 -0.111 -0.215
## Residual          8.659992              
## 
## Fixed effects: CHANGE ~ TRIAL * scale(SMR_contr, scale = FALSE) 
##                                                      Value Std.Error  DF
## (Intercept)                                       50.38244  2.385594 116
## TRIALHypoxia                                     -50.22621  2.492663 116
## TRIALLowSalinity                                 -42.22318  4.219647 116
## scale(SMR_contr, scale = FALSE)                  -21.55293  3.820913  58
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)      13.51107  3.992402 116
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)  10.18925  6.758444 116
##                                                     t-value p-value
## (Intercept)                                       21.119453  0.0000
## TRIALHypoxia                                     -20.149618  0.0000
## TRIALLowSalinity                                 -10.006330  0.0000
## scale(SMR_contr, scale = FALSE)                   -5.640779  0.0000
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)       3.384195  0.0010
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)   1.507633  0.1344
##  Correlation: 
##                                                  (Intr) TRIALH TRIALL s(Ss=F
## TRIALHypoxia                                     -0.620                     
## TRIALLowSalinity                                 -0.215 -0.035              
## scale(SMR_contr, scale = FALSE)                   0.000  0.000  0.000       
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)      0.000  0.000  0.000 -0.620
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)  0.000  0.000  0.000 -0.215
##                                                  TRIALHs=F
## TRIALHypoxia                                              
## TRIALLowSalinity                                          
## scale(SMR_contr, scale = FALSE)                           
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)              
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -0.035   
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.66585731 -0.30692682  0.02265975  0.29410064  1.72690964 
## 
## Number of Observations: 180
## Number of Groups: 60

Conclusions:

• the average change in SMR for the High Temperature trial at the average control SMR is 50.38.
• on average, the change in SMR during Hypoxia is 50.23 less at the average control SMR than during the High Temperature trial.
• on average, the change in SMR during Low Salinity is 42.22 less at the average control SMR than during the High Temperature trial.
• during the High Temperature trial, for every one unit change in control SMR, the change in SMR declined by 21.55
• the rate of change in SMR during the Hypoxia trial was 13.51 more than that associated with High Temperature.
• the rate of change in SMR during the Low Salinity trial was 10.19 more than that associated with High Temperature (although this was not found to be significant).
• there is approximately twice as much variation within Fish (replicate level) as there is between Fish (22.09 vs
  10).

confint

intervals(norin.lme3a)

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                                                        lower      est.
## (Intercept)                                       44.1837220  50.38244
## TRIALHypoxia                                     -58.2125517 -50.22621
## TRIALLowSalinity                                 -50.2095159 -42.22318
## scale(SMR_contr, scale = FALSE)                  -31.5868953 -21.55293
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)       0.7196578  13.51107
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)  -2.6021570  10.18925
##                                                      upper
## (Intercept)                                       56.58115
## TRIALHypoxia                                     -42.23987
## TRIALLowSalinity                                 -34.23684
## scale(SMR_contr, scale = FALSE)                  -11.51896
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)      26.30248
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)  22.98066
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: FISHID 
##                    lower     est.    upper
## sd((Intercept)) 5.931668 9.996282 16.84613
## 
##  Within-group standard error:
##    lower     est.    upper 
## 19.41879 22.08543 25.11827

tidy

tidy(norin.lme3a,  conf.int=TRUE)

tidy(norin.lme3a,  conf.int=TRUE) %>% kable

effect	group	term	estimate	std.error	df	statistic	p.value	conf.low	conf.high
fixed	fixed	(Intercept)	50.382438	3.129675	116	16.098295	0.0000000	44.1837220	56.58115
fixed	fixed	TRIALHypoxia	-50.226212	4.032230	116	-12.456187	0.0000000	-58.2125517	-42.23987
fixed	fixed	TRIALLowSalinity	-42.223176	4.032230	116	-10.471420	0.0000000	-50.2095159	-34.23684
fixed	fixed	scale(SMR_contr, scale = FALSE)	-21.552927	5.012679	58	-4.299682	0.0000663	-31.5868953	-11.51896
fixed	fixed	TRIALHypoxia:scale(SMR_contr, scale = FALSE)	13.511068	6.458266	116	2.092058	0.0386147	0.7196578	26.30248
fixed	fixed	TRIALLowSalinity:scale(SMR_contr, scale = FALSE)	10.189253	6.458266	116	1.577707	0.1173561	-2.6021570	22.98066
ran_pars	FISHID	sd_(Intercept)	9.996282	NA	NA	NA	NA	5.9316684	16.84613
ran_pars	Residual	sd_Observation	22.085434	NA	NA	NA	NA	NA	NA

Conclusions:

• the average change in SMR for the High Temperature trial at the average control SMR is 50.38.
• on average, the change in SMR during Hypoxia is 50.23 less at the average control SMR than during the High Temperature trial.
• on average, the change in SMR during Low Salinity is 42.22 less at the average control SMR than during the High Temperature trial.
• during the High Temperature trial, for every one unit change in control SMR, the change in SMR declined by 21.55
• the rate of change in SMR during the Hypoxia trial was 13.51 more than that associated with High Temperature.
• the rate of change in SMR during the Low Salinity trial was 10.19 more than that associated with High Temperature (although this was not found to be significant).
• there is approximately twice as much variation within Fish (replicate level) as there is between Fish (22.09 vs
  10).

tab_model

# warning this is only appropriate for html output
sjPlot::tab_model(norin.lme3a, show.se=TRUE, show.aic=TRUE)

	CHANGE
Predictors	Estimates	std. Error	CI	p
(Intercept)	50.38	3.13	44.18 – 56.58	<0.001
TRIAL [Hypoxia]	-50.23	4.03	-58.21 – -42.24	<0.001
TRIAL [LowSalinity]	-42.22	4.03	-50.21 – -34.24	<0.001
SMR_contr	-21.55	5.01	-31.59 – -11.52	<0.001
TRIAL [Hypoxia] * SMR_contr	13.51	6.46	0.72 – 26.30	0.039
TRIAL [LowSalinity] * SMR_contr	10.19	6.46	-2.60 – 22.98	0.117
Random Effects
σ2	487.77
τ00 FISHID	99.93
N FISHID	60
Observations	180
Marginal R2 / Conditional R2	0.541 / NA
AIC	1636.349

lmer (lme4)

summary

summary(norin.lmer3a)

## Linear mixed model fit by REML ['lmerMod']
## Formula: CHANGE ~ TRIAL * scale(SMR_contr, scale = FALSE) + (1 | FISHID)
##    Data: norin
## 
## REML criterion at convergence: 1620.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8551 -0.6207  0.0224  0.5421  3.7421 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  FISHID   (Intercept)  99.93    9.996  
##  Residual             487.77   22.085  
## Number of obs: 180, groups:  FISHID, 60
## 
## Fixed effects:
##                                                  Estimate Std. Error t value
## (Intercept)                                        50.382      3.130  16.098
## TRIALHypoxia                                      -50.226      4.032 -12.456
## TRIALLowSalinity                                  -42.223      4.032 -10.471
## scale(SMR_contr, scale = FALSE)                   -21.553      5.013  -4.300
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)       13.511      6.458   2.092
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)   10.189      6.458   1.578
## 
## Correlation of Fixed Effects:
##             (Intr) TRIALH TRIALL s(Ss=F TRIALHs=F
## TRIALHypoxi -0.644                               
## TRIALLwSlnt -0.644  0.500                        
## s(SMR_,s=FA  0.000  0.000  0.000                 
## TRIALH:(s=F  0.000  0.000  0.000 -0.644          
## TRIALLS:s=F  0.000  0.000  0.000 -0.644  0.500

Conclusions:

• the average change in SMR for the High Temperature trial at the average control SMR is 50.38.
• on average, the change in SMR during Hypoxia is 50.23 less at the average control SMR than during the High Temperature trial.
• on average, the change in SMR during Low Salinity is 42.22 less at the average control SMR than during the High Temperature trial.
• during the High Temperature trial, for every one unit change in control SMR, the change in SMR declined by 21.55
• the rate of change in SMR during the Hypoxia trial was 13.51 more than that associated with High Temperature.
• the rate of change in SMR during the Low Salinity trial was 10.19 more than that associated with High Temperature (although this was not found to be significant).
• there is approximately twice as much variation within Fish (replicate level) as there is between Fish (22.09 vs
  10).

confint

confint(norin.lmer3a)

##                                                        2.5 %    97.5 %
## .sig01                                             3.1895683  14.79565
## .sigma                                            19.2319058  24.78039
## (Intercept)                                       44.3172310  56.44765
## TRIALHypoxia                                     -58.0590014 -42.39342
## TRIALLowSalinity                                 -50.0559656 -34.39039
## scale(SMR_contr, scale = FALSE)                  -31.2673343 -11.83852
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)       0.9655934  26.05654
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)  -2.3562215  22.73473

tidy

tidy(norin.lmer3a,  conf.int=TRUE)

tidy(norin.lmer3a,  conf.int=TRUE) %>% kable

effect	group	term	estimate	std.error	statistic	conf.low	conf.high
fixed	NA	(Intercept)	50.382438	3.129675	16.098295	44.2483875	56.51649
fixed	NA	TRIALHypoxia	-50.226212	4.032230	-12.456186	-58.1292378	-42.32319
fixed	NA	TRIALLowSalinity	-42.223176	4.032230	-10.471420	-50.1262020	-34.32015
fixed	NA	scale(SMR_contr, scale = FALSE)	-21.552927	5.012679	-4.299682	-31.3775983	-11.72826
fixed	NA	TRIALHypoxia:scale(SMR_contr, scale = FALSE)	13.511068	6.458267	2.092058	0.8530984	26.16904
fixed	NA	TRIALLowSalinity:scale(SMR_contr, scale = FALSE)	10.189253	6.458267	1.577707	-2.4687165	22.84722
ran_pars	FISHID	sd__(Intercept)	9.996281	NA	NA	NA	NA
ran_pars	Residual	sd__Observation	22.085435	NA	NA	NA	NA

Conclusions:

• the average change in SMR for the High Temperature trial at the average control SMR is 50.38.
• on average, the change in SMR during Hypoxia is 50.23 less at the average control SMR than during the High Temperature trial.
• on average, the change in SMR during Low Salinity is 42.22 less at the average control SMR than during the High Temperature trial.
• during the High Temperature trial, for every one unit change in control SMR, the change in SMR declined by 21.55
• the rate of change in SMR during the Hypoxia trial was 13.51 more than that associated with High Temperature.
• the rate of change in SMR during the Low Salinity trial was 10.19 more than that associated with High Temperature (although this was not found to be significant).
• there is approximately twice as much variation within Fish (replicate level) as there is between Fish (22.09 vs
  10).

tab_model

# warning this is only appropriate for html output
sjPlot::tab_model(norin.lmer3a, show.se=TRUE, show.aic=TRUE)

	CHANGE
Predictors	Estimates	std. Error	CI	p
(Intercept)	50.38	3.13	44.25 – 56.52	<0.001
TRIAL [Hypoxia]	-50.23	4.03	-58.13 – -42.32	<0.001
TRIAL [LowSalinity]	-42.22	4.03	-50.13 – -34.32	<0.001
SMR_contr	-21.55	5.01	-31.38 – -11.73	<0.001
TRIAL [Hypoxia] * SMR_contr	13.51	6.46	0.85 – 26.17	0.036
TRIAL [LowSalinity] * SMR_contr	10.19	6.46	-2.47 – 22.85	0.115
Random Effects
σ2	487.77
τ00 FISHID	99.93
ICC	0.17
N FISHID	60
Observations	180
Marginal R2 / Conditional R2	0.494 / 0.580
AIC	1636.349

glmmTMB (glmmTMB)

summary

summary(norin.glmmTMB3a)

## Linear mixed model fit by REML ['lmerMod']
## Formula: CHANGE ~ TRIAL * scale(SMR_contr, scale = FALSE) + (1 | FISHID)
##    Data: norin
## 
## REML criterion at convergence: 1620.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8551 -0.6207  0.0224  0.5421  3.7421 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  FISHID   (Intercept)  99.93    9.996  
##  Residual             487.77   22.085  
## Number of obs: 180, groups:  FISHID, 60
## 
## Fixed effects:
##                                                  Estimate Std. Error t value
## (Intercept)                                        50.382      3.130  16.098
## TRIALHypoxia                                      -50.226      4.032 -12.456
## TRIALLowSalinity                                  -42.223      4.032 -10.471
## scale(SMR_contr, scale = FALSE)                   -21.553      5.013  -4.300
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)       13.511      6.458   2.092
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)   10.189      6.458   1.578
## 
## Correlation of Fixed Effects:
##             (Intr) TRIALH TRIALL s(Ss=F TRIALHs=F
## TRIALHypoxi -0.644                               
## TRIALLwSlnt -0.644  0.500                        
## s(SMR_,s=FA  0.000  0.000  0.000                 
## TRIALH:(s=F  0.000  0.000  0.000 -0.644          
## TRIALLS:s=F  0.000  0.000  0.000 -0.644  0.500

Conclusions:

• the average change in SMR for the High Temperature trial at the average control SMR is 50.38.
• on average, the change in SMR during Hypoxia is 50.23 less at the average control SMR than during the High Temperature trial.
• on average, the change in SMR during Low Salinity is 42.22 less at the average control SMR than during the High Temperature trial.
• during the High Temperature trial, for every one unit change in control SMR, the change in SMR declined by 21.55
• the rate of change in SMR during the Hypoxia trial was 13.51 more than that associated with High Temperature.
• the rate of change in SMR during the Low Salinity trial was 10.19 more than that associated with High Temperature (although this was not found to be significant).
• there is approximately twice as much variation within Fish (replicate level) as there is between Fish (22.09 vs
  10).

confint

confint(norin.glmmTMB3a)

##                                                        2.5 %    97.5 %
## .sig01                                             3.1895683  14.79565
## .sigma                                            19.2319058  24.78039
## (Intercept)                                       44.3172310  56.44765
## TRIALHypoxia                                     -58.0590014 -42.39342
## TRIALLowSalinity                                 -50.0559656 -34.39039
## scale(SMR_contr, scale = FALSE)                  -31.2673343 -11.83852
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)       0.9655934  26.05654
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)  -2.3562215  22.73473

tidy

tidy(norin.glmmTMB3a,  conf.int=TRUE)

tidy(norin.glmmTMB3a,  conf.int=TRUE) %>% kable

effect	group	term	estimate	std.error	statistic	conf.low	conf.high
fixed	NA	(Intercept)	50.382438	3.129675	16.098295	44.2483875	56.51649
fixed	NA	TRIALHypoxia	-50.226212	4.032230	-12.456186	-58.1292378	-42.32319
fixed	NA	TRIALLowSalinity	-42.223176	4.032230	-10.471420	-50.1262020	-34.32015
fixed	NA	scale(SMR_contr, scale = FALSE)	-21.552927	5.012679	-4.299682	-31.3775983	-11.72826
fixed	NA	TRIALHypoxia:scale(SMR_contr, scale = FALSE)	13.511068	6.458267	2.092058	0.8530984	26.16904
fixed	NA	TRIALLowSalinity:scale(SMR_contr, scale = FALSE)	10.189253	6.458267	1.577707	-2.4687165	22.84722
ran_pars	FISHID	sd__(Intercept)	9.996281	NA	NA	NA	NA
ran_pars	Residual	sd__Observation	22.085435	NA	NA	NA	NA

Conclusions:

• the average change in SMR for the High Temperature trial at the average control SMR is 50.38.
• on average, the change in SMR during Hypoxia is 50.23 less at the average control SMR than during the High Temperature trial.
• on average, the change in SMR during Low Salinity is 42.22 less at the average control SMR than during the High Temperature trial.
• during the High Temperature trial, for every one unit change in control SMR, the change in SMR declined by 21.55
• the rate of change in SMR during the Hypoxia trial was 13.51 more than that associated with High Temperature.
• the rate of change in SMR during the Low Salinity trial was 10.19 more than that associated with High Temperature (although this was not found to be significant).
• there is approximately twice as much variation within Fish (replicate level) as there is between Fish (22.09 vs
  10).

tab_model

# warning this is only appropriate for html output
sjPlot::tab_model(norin.glmmTMB3a, show.se=TRUE, show.aic=TRUE)

	CHANGE
Predictors	Estimates	std. Error	CI	p
(Intercept)	50.38	3.13	44.25 – 56.52	<0.001
TRIAL [Hypoxia]	-50.23	4.03	-58.13 – -42.32	<0.001
TRIAL [LowSalinity]	-42.22	4.03	-50.13 – -34.32	<0.001
SMR_contr	-21.55	5.01	-31.38 – -11.73	<0.001
TRIAL [Hypoxia] * SMR_contr	13.51	6.46	0.85 – 26.17	0.036
TRIAL [LowSalinity] * SMR_contr	10.19	6.46	-2.47 – 22.85	0.115
Random Effects
σ2	487.77
τ00 FISHID	99.93
ICC	0.17
N FISHID	60
Observations	180
Marginal R2 / Conditional R2	0.494 / 0.580
AIC	1636.349

Predictions / further analyses

The presence of a significant interaction suggests that the nature of the relationship between SMR change and the control SMR differs between the three trials (High Temperature, Hypoxia and Low Salinity). Equivalently, the diffences between the three trials is not constant over all levels of control SMR.

To further tease these interactions apart, we could either:

• compare the slopes separately for each of the three trials
• compare the effects of the three trials at multiple values of control SMR

We will now explore each of the above, although typically we would only do one or the other (depending on which variable was of more interest).

lme (nlme)

Compare slopes

emtrends(norin.lme3a, pairwise~TRIAL, var='SMR_contr')

## $emtrends
##  TRIAL           SMR_contr.trend   SE df lower.CL upper.CL
##  HighTemperature          -21.55 5.01 58    -31.6   -11.52
##  Hypoxia                   -8.04 5.01 58    -18.1     1.99
##  LowSalinity              -11.36 5.01 58    -21.4    -1.33
## 
## Degrees-of-freedom method: containment 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                      estimate   SE  df t.ratio p.value
##  HighTemperature - Hypoxia       -13.51 6.46 116 -2.092  0.0959 
##  HighTemperature - LowSalinity   -10.19 6.46 116 -1.578  0.2594 
##  Hypoxia - LowSalinity             3.32 6.46 116  0.514  0.8645 
## 
## Degrees-of-freedom method: containment 
## P value adjustment: tukey method for comparing a family of 3 estimates

Conclusions:

• the partial slope associated with control SMR during High Temperature was found to be 13.51 units more negative than during Hypoxia (although not significant). Note, we already had this partial slope in the original summary table.
• the partial slope associated with control SMR during High Temperature was found to be 10.19 units more negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.
• the partial slope associated with control SMR during Hypoxia was found to be 3.32 units less negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.

Compare Trials (at different values of control SMR)

We will estimated the pairwise differences between each of the three trials at the minimum, mean and maximum control SMR.

norin.grid <- with(norin,  list(SMR_contr=c(min(SMR_contr), mean(SMR_contr), max(SMR_contr))))
norin.grid

## $SMR_contr
## [1] 3.984319 5.178857 6.700537

emmeans(norin.lme3a, pairwise~TRIAL|SMR_contr,  at=norin.grid)

## $emmeans
## SMR_contr = 3.98:
##  TRIAL            emmean   SE df lower.CL upper.CL
##  HighTemperature  76.128 6.76 58    62.60    89.65
##  Hypoxia           9.763 6.76 58    -3.76    23.29
##  LowSalinity      21.734 6.76 58     8.21    35.26
## 
## SMR_contr = 5.18:
##  TRIAL            emmean   SE df lower.CL upper.CL
##  HighTemperature  50.382 3.13 59    44.12    56.64
##  Hypoxia           0.156 3.13 59    -6.11     6.42
##  LowSalinity       8.159 3.13 59     1.90    14.42
## 
## SMR_contr = 6.70:
##  TRIAL            emmean   SE df lower.CL upper.CL
##  HighTemperature  17.586 8.24 58     1.08    34.09
##  Hypoxia         -12.081 8.24 58   -28.58     4.42
##  LowSalinity      -9.133 8.24 58   -25.64     7.37
## 
## Degrees-of-freedom method: containment 
## Confidence level used: 0.95 
## 
## $contrasts
## SMR_contr = 3.98:
##  contrast                      estimate    SE  df t.ratio p.value
##  HighTemperature - Hypoxia        66.37  8.70 116  7.624  <.0001 
##  HighTemperature - LowSalinity    54.39  8.70 116  6.249  <.0001 
##  Hypoxia - LowSalinity           -11.97  8.70 116 -1.375  0.3572 
## 
## SMR_contr = 5.18:
##  contrast                      estimate    SE  df t.ratio p.value
##  HighTemperature - Hypoxia        50.23  4.03 116 12.456  <.0001 
##  HighTemperature - LowSalinity    42.22  4.03 116 10.471  <.0001 
##  Hypoxia - LowSalinity            -8.00  4.03 116 -1.985  0.1205 
## 
## SMR_contr = 6.70:
##  contrast                      estimate    SE  df t.ratio p.value
##  HighTemperature - Hypoxia        29.67 10.62 116  2.793  0.0167 
##  HighTemperature - LowSalinity    26.72 10.62 116  2.515  0.0351 
##  Hypoxia - LowSalinity            -2.95 10.62 116 -0.278  0.9584 
## 
## Degrees-of-freedom method: containment 
## P value adjustment: tukey method for comparing a family of 3 estimates

Conclusions:

• as control SMR increases, the magnitude of differences between the trials becomes reduced.
• the High Temperature trial resulted in the greatest positive change in SMR.
• the change in SMR was not found to be different between the Hypoxia and Low Salinity trials.

\(R^2\)

r.squaredGLMM(norin.lme3a)

##            R2m       R2c
## [1,] 0.4942091 0.5802091

## Nakagawa's R2
performance::r2_nakagawa(norin.lme3a)

## Random effect variances not available. Returned R2 does not account for random effects.

## # R2 for Mixed Models
## 
##   Conditional R2: NA
##      Marginal R2: 0.541

lmer (lmer)

Compare slopes

emtrends(norin.lmer3a, pairwise~TRIAL, var='SMR_contr')

## $emtrends
##  TRIAL           SMR_contr.trend   SE  df lower.CL upper.CL
##  HighTemperature          -21.55 5.01 164    -31.5   -11.66
##  Hypoxia                   -8.04 5.01 164    -17.9     1.86
##  LowSalinity              -11.36 5.01 164    -21.3    -1.47
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                      estimate   SE  df t.ratio p.value
##  HighTemperature - Hypoxia       -13.51 6.46 116 -2.092  0.0959 
##  HighTemperature - LowSalinity   -10.19 6.46 116 -1.578  0.2594 
##  Hypoxia - LowSalinity             3.32 6.46 116  0.514  0.8645 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

norin.emt <- emtrends(norin.lmer3a, pairwise~TRIAL, var='SMR_contr')$contrasts %>% as.data.frame

Conclusions:

• the partial slope associated with control SMR during High Temperature was found to be 13.51 units more negative than during Hypoxia (although not significant). Note, we already had this partial slope in the original summary table.
• the partial slope associated with control SMR during High Temperature was found to be 10.19 units more negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.
• the partial slope associated with control SMR during Hypoxia was found to be 3.32 units less negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.

Compare Trials (at different values of control SMR)

norin.grid <- with(norin,  list(SMR_contr=c(min(SMR_contr), mean(SMR_contr), max(SMR_contr))))
norin.grid

## $SMR_contr
## [1] 3.984319 5.178857 6.700537

emmeans(norin.lmer3a, pairwise~TRIAL|SMR_contr,  at=norin.grid)

## $emmeans
## SMR_contr = 3.98:
##  TRIAL            emmean   SE  df lower.CL upper.CL
##  HighTemperature  76.128 6.76 164    62.79    89.47
##  Hypoxia           9.763 6.76 164    -3.58    23.10
##  LowSalinity      21.734 6.76 164     8.39    35.07
## 
## SMR_contr = 5.18:
##  TRIAL            emmean   SE  df lower.CL upper.CL
##  HighTemperature  50.382 3.13 164    44.20    56.56
##  Hypoxia           0.156 3.13 164    -6.02     6.34
##  LowSalinity       8.159 3.13 164     1.98    14.34
## 
## SMR_contr = 6.70:
##  TRIAL            emmean   SE  df lower.CL upper.CL
##  HighTemperature  17.586 8.24 164     1.31    33.87
##  Hypoxia         -12.081 8.24 164   -28.36     4.20
##  LowSalinity      -9.133 8.24 164   -25.41     7.15
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
## SMR_contr = 3.98:
##  contrast                      estimate    SE  df t.ratio p.value
##  HighTemperature - Hypoxia        66.37  8.70 116  7.624  <.0001 
##  HighTemperature - LowSalinity    54.39  8.70 116  6.249  <.0001 
##  Hypoxia - LowSalinity           -11.97  8.70 116 -1.375  0.3572 
## 
## SMR_contr = 5.18:
##  contrast                      estimate    SE  df t.ratio p.value
##  HighTemperature - Hypoxia        50.23  4.03 116 12.456  <.0001 
##  HighTemperature - LowSalinity    42.22  4.03 116 10.471  <.0001 
##  Hypoxia - LowSalinity            -8.00  4.03 116 -1.985  0.1205 
## 
## SMR_contr = 6.70:
##  contrast                      estimate    SE  df t.ratio p.value
##  HighTemperature - Hypoxia        29.67 10.62 116  2.793  0.0167 
##  HighTemperature - LowSalinity    26.72 10.62 116  2.515  0.0351 
##  Hypoxia - LowSalinity            -2.95 10.62 116 -0.278  0.9584 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

Conclusions:

• as control SMR increases, the magnitude of differences between the trials becomes reduced.
• the High Temperature trial resulted in the greatest positive change in SMR.
• the change in SMR was not found to be different between the Hypoxia and Low Salinity trials.

\(R^2\)

r.squaredGLMM(norin.lmer3a)

##            R2m      R2c
## [1,] 0.4942091 0.580209

## Nakagawa's R2
performance::r2_nakagawa(norin.lmer3a)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.580
##      Marginal R2: 0.494

glmmTMB (glmmTMB)

Compare slopes

emtrends(norin.glmmTMB3a, pairwise~TRIAL, var='SMR_contr')

## $emtrends
##  TRIAL           SMR_contr.trend   SE  df lower.CL upper.CL
##  HighTemperature          -21.55 5.01 164    -31.5   -11.66
##  Hypoxia                   -8.04 5.01 164    -17.9     1.86
##  LowSalinity              -11.36 5.01 164    -21.3    -1.47
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                      estimate   SE  df t.ratio p.value
##  HighTemperature - Hypoxia       -13.51 6.46 116 -2.092  0.0959 
##  HighTemperature - LowSalinity   -10.19 6.46 116 -1.578  0.2594 
##  Hypoxia - LowSalinity             3.32 6.46 116  0.514  0.8645 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

norin.emt <- emtrends(norin.glmmTMB3a, pairwise~TRIAL, var='SMR_contr')$contrasts %>% as.data.frame

Conclusions:

• the partial slope associated with control SMR during High Temperature was found to be 13.51 units more negative than during Hypoxia (although not significant). Note, we already had this partial slope in the original summary table.
• the partial slope associated with control SMR during High Temperature was found to be 10.19 units more negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.
• the partial slope associated with control SMR during Hypoxia was found to be 3.32 units less negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.

Compare Trials (at different values of control SMR)

norin.grid <- with(norin,  list(SMR_contr=c(min(SMR_contr), mean(SMR_contr), max(SMR_contr))))
norin.grid

## $SMR_contr
## [1] 3.984319 5.178857 6.700537

emmeans(norin.glmmTMB3a, pairwise~TRIAL|SMR_contr,  at=norin.grid)

## $emmeans
## SMR_contr = 3.98:
##  TRIAL            emmean   SE  df lower.CL upper.CL
##  HighTemperature  76.128 6.76 164    62.79    89.47
##  Hypoxia           9.763 6.76 164    -3.58    23.10
##  LowSalinity      21.734 6.76 164     8.39    35.07
## 
## SMR_contr = 5.18:
##  TRIAL            emmean   SE  df lower.CL upper.CL
##  HighTemperature  50.382 3.13 164    44.20    56.56
##  Hypoxia           0.156 3.13 164    -6.02     6.34
##  LowSalinity       8.159 3.13 164     1.98    14.34
## 
## SMR_contr = 6.70:
##  TRIAL            emmean   SE  df lower.CL upper.CL
##  HighTemperature  17.586 8.24 164     1.31    33.87
##  Hypoxia         -12.081 8.24 164   -28.36     4.20
##  LowSalinity      -9.133 8.24 164   -25.41     7.15
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
## SMR_contr = 3.98:
##  contrast                      estimate    SE  df t.ratio p.value
##  HighTemperature - Hypoxia        66.37  8.70 116  7.624  <.0001 
##  HighTemperature - LowSalinity    54.39  8.70 116  6.249  <.0001 
##  Hypoxia - LowSalinity           -11.97  8.70 116 -1.375  0.3572 
## 
## SMR_contr = 5.18:
##  contrast                      estimate    SE  df t.ratio p.value
##  HighTemperature - Hypoxia        50.23  4.03 116 12.456  <.0001 
##  HighTemperature - LowSalinity    42.22  4.03 116 10.471  <.0001 
##  Hypoxia - LowSalinity            -8.00  4.03 116 -1.985  0.1205 
## 
## SMR_contr = 6.70:
##  contrast                      estimate    SE  df t.ratio p.value
##  HighTemperature - Hypoxia        29.67 10.62 116  2.793  0.0167 
##  HighTemperature - LowSalinity    26.72 10.62 116  2.515  0.0351 
##  Hypoxia - LowSalinity            -2.95 10.62 116 -0.278  0.9584 
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

Conclusions:

• as control SMR increases, the magnitude of differences between the trials becomes reduced.
• the High Temperature trial resulted in the greatest positive change in SMR.
• the change in SMR was not found to be different between the Hypoxia and Low Salinity trials.

\(R^2\)

r.squaredGLMM(norin.glmmTMB3a)

##            R2m      R2c
## [1,] 0.4942091 0.580209

## Nakagawa's R2
performance::r2_nakagawa(norin.glmmTMB3a)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.580
##      Marginal R2: 0.494

Summary figures

lme (nlme)

norin.grid=with(norin, list(SMR_contr=seq(min(SMR_contr), max(SMR_contr), len=100)))
newdata = emmeans(norin.lme3a, ~SMR_contr|TRIAL, at=norin.grid) %>%
    as.data.frame
head(newdata)

ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
  geom_ribbon(aes(ymin=lower.CL, ymax=upper.CL, fill=TRIAL), alpha=0.3) +
  geom_line(aes(, color=TRIAL)) +
  theme_classic() +
  theme(legend.position = c(0.99, 0.99),
        legend.justification = c(1, 1))

obs <- norin.lme3a %>%
  augment() %>%
  mutate(PartialObs=.fitted + .resid)

ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
  geom_ribbon(aes(ymin=lower.CL, ymax=upper.CL, fill=TRIAL), alpha=0.3) +
  geom_line(aes(, color=TRIAL)) +
  geom_point(data=obs,  aes(y=PartialObs,  color=TRIAL)) +
  theme_classic()

lmer (lme4)

norin.grid=with(norin, list(SMR_contr=seq(min(SMR_contr), max(SMR_contr), len=100)))
newdata = emmeans(norin.lmer3a, ~SMR_contr|TRIAL, at=norin.grid) %>%
    as.data.frame
head(newdata)

ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
  geom_ribbon(aes(ymin=lower.CL, ymax=upper.CL, fill=TRIAL), alpha=0.3) +
  geom_line(aes(, color=TRIAL)) +
  theme_classic() +
  theme(legend.position = c(0.99, 0.99),
        legend.justification = c(1, 1))

obs <- norin.lmer3a %>%
  augment() %>%
  bind_cols(norin %>% select(SMR_contr)) %>%
  mutate(PartialObs=.fitted + .resid)

ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
  geom_ribbon(aes(ymin=lower.CL, ymax=upper.CL, fill=TRIAL), alpha=0.3) +
  geom_line(aes(, color=TRIAL)) +
  geom_point(data=obs,  aes(y=PartialObs,  color=TRIAL)) +
  theme_classic()

glmmTMB (glmmTMB)

norin.grid=with(norin, list(SMR_contr=seq(min(SMR_contr), max(SMR_contr), len=100)))
newdata = emmeans(norin.glmmTMB3a, ~SMR_contr|TRIAL, at=norin.grid) %>%
    as.data.frame
head(newdata)

ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
  geom_ribbon(aes(ymin=lower.CL, ymax=upper.CL, fill=TRIAL), alpha=0.3) +
  geom_line(aes(, color=TRIAL)) +
  theme_classic() +
  theme(legend.position = c(0.99, 0.99),
        legend.justification = c(1, 1))

obs <- norin.glmmTMB3a %>%
  augment() %>%
  bind_cols(norin %>% select(SMR_contr)) %>%
  mutate(PartialObs=.fitted + .resid)

ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
  geom_ribbon(aes(ymin=lower.CL, ymax=upper.CL, fill=TRIAL), alpha=0.3) +
  geom_line(aes(, color=TRIAL)) +
  geom_point(data=obs,  aes(y=PartialObs,  color=TRIAL)) +
  theme_classic()

References

Norin, Tommy, Hans Malte, and Timothy D. Clark. 2015. “Differential Plasticity of Metabolic Rate Phenotypes in a Tropical Fish Facing Environmental Change.” Functional Ecology 30 (3): 369–78. https://doi.org/10.1111/1365-2435.12503.